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Abstract 

We study the interaction between the Be star and the pulsar in the TeV binary PSR B1259— 63/LS 2883, 
using 3-D SPH simulations of the tidal and wind interactions in this Be-pulsar system. We first run a 
simulation without pulsar wind nor Be wind, taking into account only the gravitational effect of the pulsar 
on the Be disk. In this simulation, the gas particles are ejected at a constant rate from the equatorial 
surface of the Be star, which is tilted in a direction consistent with multi-waveband observations. We run 
the simulation until the Be disk is fully developed and starts to repeat a regular tidal interaction with the 
pulsar. Then, we turn on the pulsar wind and the Be wind. We run two simulations with different wind 
mass-loss rates for the Be star, one for a B2V type and the other for a significantly earlier spectral type. 
Although the global shape of the interaction surface between the pulsar wind and the Be wind agrees with 
the analytical solution, the effect of the pulsar wind on the Be disk is profound. The pulsar wind strips off 
an outer part of the Be disk, truncating the disk at a radius significantly smaller than the pulsar orbit. Our 
results, therefore, rule out the idea that the pulsar passes through the Be disk around periastron, which 
has been assumed in the previous studies. It also turns out that the location of the contact discontinuity 
can be significantly different between phases when the pulsar wind directly hits the Be disk and those 
when the pulsar wind collides with the Be wind. It is thus important to adequately take into account the 
circumstellar environment of the Be star, in order to construct a satisfactory model for this prototypical 
TeV binary. 

Keywords: gamma rays: theory - stars: emission-line, Be - stars: winds, outflows - stars: individual 
(B 1259-63) 



1. Introduction 

PSR B1259— 63/LS 2883 is a massive binary system 
consisting of a 48-ms radio pulsar and a Be star with a cir- 
cumstellar disk (Johnston et al. 1992). The system is one 
of only three binaries from which periodic TeV gamma-ray 
emission has been detected (the others are LS I +61—303 
and LS 5039; e.g., Paredes 2008). Among these TeV 
gamma-ray binaries, PSR B1259— 63/LS 2883 has several 
distinctive properties. First, the compact star of this sys- 
tem is identified as a radio pulsar, namely neutron star, 
while the nature of the compact object of the other two 
sources is still under debate. Second, the orbital period 



(3.4 yr) is very long compared to the other systems (3.9 d 
for LS 5039 and 26.5 d for LS I +61 303), and the eccen- 
tricity is so high that the binary separation at periastron is 
less than a tenth of that at apastron (Johnston et al. 1994). 
Third, nonthermal emission is detected in X-ray and radio 
wave bands as well, and light curves of all bands includ- 
ing TeV gamma-rays exhibit multiple-peaked features, in 
most cases one before periastron and another after it, al- 
though the curves differ orbit to orbit (Chernyakova et al. 
2009; Aharonian et al. 2009). 

Considering the passage of the neutron star through a 
misaligned Be disk, shock acceleration in the disk-pulsar 
interaction region has been proposed as the mechanism 
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for double-peaked light curves. Modeling of gamma-ray 
emission was first discussed by Tavani et al. (1994) for 
McV energies following the CGRO observation. They 
considered shock acceleration in the interaction region, 
and showed that the double-peaked feature appeared in 
X-ray and MeV gamma-ray ranges. Much theoretical 
work for high-energy emissions followed (e.g. Tavani & 
Arons 1997; Ball k Kirk 2000; Murata et al. 2003) and 
the emission in TeV energy range has been discussed in 
Ball &; Kirk (2000) in particular. Assuming a disk-like 
outflow around the Be star, Kawachi et al. (2004) pro- 
posed a double-peaked light curve in TeV range, which 
was followed by the TeV gamma-ray detection with a 
marginal pre-periastron peak and a clear post-periastron 
peak (Aharonian et al. 2005). As shown in the radio ob- 
servations, the light curve varies from orbit to orbit and 
the TeV double-peak was not clearly seen in the next 2007 
pcriastron passage. Aharonian et al. (2009) suggests that 
the detection of the TeV emission 50 days prior to the 
periastron disfavors that the pulsar-disk interaction is the 
primary TeV emission mechanism. 

It can be pointed out that all the previous theoretical 
studies have either neglected the recent progress in Be star 
research and adopted an out-of-dated view of the Be disk 
or neglected even the presence of the Be disk. Given that 
the gas pressure in the Be disk is much higher than the 
ram pressure of the Be wind, it is important to take into 
account the interaction between the pulsar wind and the 
Be disk adequately In this series of papers, we will study 
the high energy emission from PSR B1259-63/LS 2883, 
for the first time based on 3-D numerical simulations with 
the latest model of the Be disk. As a first step, this paper 
explores the hydrodynamic interaction between the pulsar 
wind and the circumstcllar environment of the Be star. 

The structure of the paper is as follows. In section 2, 
we first describe the Be disk model and the numerical 
method. We then discuss the viscous evolution of the Be 
disk under the tidal interaction with the pulsar, based on 
simulations where only the tidal interaction is taken into 
account. In section 3, we show the results from simula- 
tions that also take account of both the pulsar wind and 
the Be wind, and examine the effects of the pulsar wind on 
the circumstellar environment of the Be star. In section 4, 
we discuss implications for multi-wavelength observations 
and summarize the conclusions. 

2. Tidal Effect of the Pulsar on the Be Disk 

2.1. Viscous decretion disk model for Be stars 

In this subsection, we briefly summarize key observa- 
tional features of Be stars and describe the up-to-date 
model of the Be disk that is widely accepted in the Be star 
community (see Porter & Rivinius 2003 for more detailed 
descriptions of Be stars and their circumstellar disks) . 

A Be star has a two-component extended atmosphere, a 
polar wind region and an equatorial disk region. The for- 
mer consists of a low-density fast flow (~ 10 3 kms _1 ) emit- 
ting UV radiation. The wind structure is well explained 
by the line-driven wind model (Castor et al. 1975; Friend 



& Abbott 1986). On the other hand, the equatorial disk 
consists of a high-density plasma, from which the optical 
emission lines and the IR excess arise. The relationship 
between the disk size resolved with the optical interferom- 
eters (Quirrenbach et al. 1994; Quirrenbach et al. 1997) 
and the separation of the two peaks of the Ha line profile 
is in agreement with that expected for a Keplerian disk. 
The outflow in the Be disk is very subsonic. Hanuschik 
(1994, 2000) and Waters & Marlborough (1994) showed 
that the radial velocity of the disk is smaller than a few 
kms -1 , at least within ~ 10 stellar radii. Thus, it is un- 
likely that the Be disk is formed by channeling/focusing 
of the polar wind. Indeed, it is established that the Be 
disk cannot be modeled as a low velocity wind. 

One model which can reproduce all key features of the 
Be disk and is thus widely accepted, is the viscous decre- 
tion disk model proposed by Lee et al. (1991) and devel- 
oped by Porter (1999) and others (Okazaki 2001, 2007; 
Carciofi & Bjorkman 2006; Carciofi 2010). The model as- 
sumes that the star can eject material with the Keplerian 
velocity at the stellar equatorial surface. The ejected ma- 
terial then drifts outward by viscous diffusion and forms 
a geometrically thin Keplerian disk. Recent Monte Carlo 
rcdiativc transfer simulations show that such a disk has 
a temperature approximately constant at ~ 60% of the 
stellar effective temperature, except in an inner optically- 
thick region where the disk is significantly cooler (Carciofi 
& Bjorkman 2006). It should be noted that the viscous 
decretion disk model also provides a firm basis for the 
study of tidal interaction between the Be disk and the 
neutron star in Be/X-ray binaries (Negueruela & Okazaki 
2001; Okazaki & Negueruela 2001). 

2.2. Numerical model 

The simulation presented below was performed with 
a three dimensional Smoothed Particle Hydrodynamics 
(SPH) code. The code is basically the same as that used 
by Okazaki et al. (2002) (see also Bate et al. 1995). Using 
a variable smoothing length, the SPH equations with a 
standard cubic-spline kernel are integrated with an indi- 
vidual time step for each particle. In our code, the Be 
disk is modeled by an ensemble of gas particles with neg- 
ligible self-gravity. For simplicity, the gas particles are 
assumed to be isothermal at 0.6 T ff with T c r being the 
effective temperature of the Be star. On the other hand, 
the Be star and the neutron star are represented by sink 
particles with the appropriate gravitational mass. Gas 
particles that fall within a specified accretion radius are 
accreted by the sink particle. 

In simulations shown in this section, the numerical vis- 
cosity is adjusted so as to keep the Shakura-Sunyaev vis- 
cosity parameter ass = 0.1, using the approximate relation 
a S s = O.laspH V# and #3PH = (Okazaki et al. 2002), 
where aspH and /?sph are the artificial viscosity parame- 
ters 1 , and h and H are the smoothing length of individual 

The artificial viscosity commonly used in SPH consists of two 
terms: a term that is linear in the velocity differences between 
particles, which produces a shear and bulk viscosity, and a term 
that is quadratic in the velocity differences, which is needed to 
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polar axis 
z 45 deg. / 



Table 1. Model parameters 




Be star parameters 



pulsar 



Fig. 1. The model configuration of the VHE gamma-ray bi- 
nary PSR B1259-63/LS 2883. The binary orbit is in the x-y 
plane. The polar axis of Be star is tilted by 45° to the £-axis 
and has 19° of azimuth of tilt to the x-axis. 

particles and the scale-height of the Be decretion disk, re- 
spectively. No effect of pulsar wind or Be wind is taken 
into account in these simulations. 

We set the binary orbit in the x-y plane with the ma- 
jor axis along the x-axis. At t = 0, the pulsar is at the 
apastron (Phase ~ 0.5). It orbits about the Be star with 
the orbital period P or b = 1236.79 d and the eccentricity 
e = 0.87. It has long been expected that the Be disk is 
inclined from the orbital plane, because the pulsed radio 
emission disappears for about five weeks centered on the 
epoch of periastron, e.g., ~ 18 d prior to the 2004 perias- 
tron through ~ 16 d after it (Johnston et al. 2005). We set 
the tilt angle between the binary orbital axis and the Be 
star's polar axis to 45° as a parameter and the azimuth 
of tilt, i.e., the azimuthal angle of the Be star's polar axis 
from the z-axis (the direction of apastron) to 19° as sug- 
gested in Chernyakova et al. (2006) (see Fig. 1). In order 
to emulate the mass ejection from a Be star, we inject 
gas particles just outside the stellar equatorial surface at 
the rate of 3.5 x lO _9 M0yr _1 , which gives rise to a typi- 
cal disk base density of 10 _11 g cm" 3 . Once injected, gas 
particles interact with each other. As a result, most in- 
jected particles fall onto the Be star by losing angular mo- 
mentum, but a small fraction of particles drift outwards, 
getting the angular momentum from the other particles, 
and form a disk. 

Table 1 summarizes the parameters adopted in the 
following simulations. Note that not all parameters in 
Table 1 are independent. For instance, all Be disk pa- 
rameters except for the tilt angles are derived from other 
parameters. 

2.3. Long-term evolution of the Be disk under no influ- 
ence of the pulsar wind 

Figure 2 shows the evolution of surface density from t = 
to 12P or b (panel a) and the disk structure at t = 12P or b 
(panel b). In Fig. 2, the volume density is integrated verti- 
cally and averaged azimuthally, while the velocity compo- 
nents are averaged vertically and azimuthally. Figure 2(a) 
exhibits how a decretion disk forms if there is no in- 
fluence of the pulsar wind. The disk density gradually 



Mass of the Be star M* = 10 M a 
Radius of the Be star = 6 P Q a 
Effective temperature T cS = 22, 800 K b 
Critical velocity V CT n(R*) = 563.9kms~ 1 
Wind velocity V^ind = 10 3 km s _1 
Wind mass-loss rate 

10~ 9 MQyr -1 ("weak" case) 
lO _8 M0yr _1 ( "strong" case) 



eliminate particle intcrpenetration in high Mach number shocks. 
The parameters osph an d /?SPH control the linear and quadratic 
terms, respectively. 



Mass injection rate to the disk 

Mdisk = 3.5 x lO _9 M yr- lc 

Be disk parameters 
Disk temperature = 13, 680 K d 
Disk thickness P(P*)/P* = 0.024 
Sound speed c s = 13.8kms _1 

Initial tilt angles: Misalignment by 45° from the bi- 
nary orbital axis, with the polar axis projection of 
19° from apastron c (see Fig. 1) 
Pulsar and orbital parameters 
Pulsar wind power Ppsr = 8.2 x 10 35 ergs _1 a 
Orbital period P orb = 1236.79d (= 3.386yr) f 
Orbital eccentricity e = 0.87 a 
Mass ratio q = 0.14 (M PS r = 1.4M0) a 
Semimajor axis a = 181.82^ = 7.59 ■ 10 13 cm a 
a Johnston et al. (1994) 

b The effective temperature of a B2V star (Cox 2000) 
c Mass ejection rate that gives the typical disk base 
density of 10~ n g cm 3 

d Assumed to be isothermal at 0.6T e ff (Carciofi & 

Bjorkman 2006). 

e Chernyakova et al. (2006) 

f Taken from Johnston et al. 1994. Recent timing 
analyses give P orb = 1236.72 d (Wang et al. 2004). 

increases with time and approaches an asymptotic dis- 
tribution. The timescale of the disk growth is the vis- 
cous timescale, which is given by [ass(PA) 2 ^K] -1 , where 
= (G-/VP/P 3 ) 1 / 2 is the Keplerian angular frequency. 
For parameters summarized in Table 1, this timescale is 
~12(a/0.1)- 1 (r/10P») 1/2 yr. Thus, at t ~ 12P orb ~ 40yr, 
the disk is already fully developed, as seen in Fig. 2(a). 

We also note from Fig. 2(a) that the disk extends be- 
yond the periastron separation (r ~ 24P*) with no break 
in the radial density distribution. This suggests that the 
Be disk in long-period, highly-eccentric systems like PSR 
B1259— 63 is not truncated solely by the tidal torques of 
the neutron star, unlike disks in other Bc-neutron star 
binaries with shorter orbital periods and lower eccentric- 
ities. In such systems, the neutron star passes through, 
and thus strongly disturbs, the Be disk around every pe- 
riastron passage, as we will see in section 2.4. 

As discussed in section 2.1, the Be disk is Keplerian. 
Fig. 2(b) shows that the angular velocity (the dashed line) 
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Fig. 2. Left: Surface density evolution from t =0 to 12 in 
units of P or b in a simulation where only the tidal interaction 
is taken into account. The interval of time between adjacent 
contours is 2 (t = 2, 4, 12 from bottom). Right: Disk 
structure at t = 12. The thick solid, dashed, and dash-dotted 
lines denote the surface density, the azimuthal velocity nor- 
malized by the stellar critical velocity, and the radial Mach 
number, respectively. For comparison, the Keplerian rotation 
distribution is shown by the thin solid line. In both panels, 
the density is integrated vertically and averaged azimuthally, 
while the velocity components arc averaged vertically and az- 
imuthally. 

for r < 100-R* is indistinguishable from the Keplerian ro- 
tation (the thin solid line). The radial velocity, V r , is very 
subsonic. It increases with radius, but even at 100i?», V r 
is still of the order of ~ 0.1 c s ~ lkms , where c s is the 
isothermal sound speed in the Be disk. These are typical 
features of viscous decretion disks around Be stars (e.g., 
Lee et al. 1991). 

2.4- Tidal interaction around periastron passage 

Figure 3 provides snapshots over a month from 13 d be- 
fore periastron passage (0 = 0.99) to 23 d after it ((f) = 
0.02). The upper and lower panels show respectively the 
column density in cgs units along the binary orbital axis 
(z-axis) and the minor axis of the binary orbit (y-axis). 
The first disk-passing event of the pulsar starts about 20 
days prior to periastron passage, while the second event 
starts ~ 12 days after it. Each event lasts for ~ 2 weeks. 
As seen in the figure, the Be disk is strongly disturbed 
by these events. The disk is temporarily warped at each 
of these events via the tidal interaction with the pulsar. 
At the same time, the tidal interaction around periastron 
excites a two-armed spiral wave, via which angular mo- 
mentum is removed from the disk (see, e.g., Artymowicz 
& Lubow 1994). This angular momentum transfer causes 
shrinking of the Be disk. Its radius gradually recovers af- 
terward by viscous diffusion (Okazaki et al. 2002). Note 
that because of the cumulative effect of the tidal torques, 
the disturbance in the Be disk is larger after periastron 
passage, despite that the distance between the Be disk 
and pulsar is closest at the first disk-passing event. 



3. Effects of the Pulsar Wind on the 
Circumstellar Environment of the Be star 

3.1. Numerical setup 

In order to see how the pulsar wind interacts with the 
Be disk and wind, we carried out simulations where both 
the pulsar and Be winds are now taken into account, us- 
ing the data from the above tidal simulation. In these new 
wind simulations, we turned on the pulsar and Be winds 
at a certain time, after the Be disk was fully developed 
in the tidal simulation. We started the wind simulation 
at t = ll.44.Port,, 74 days prior to periastron, much longer 
than the crossing time of the Be wind over the simula- 
tion volume, ~ 10 d. The relativistic pulsar wind crosses 
the simulation volume much faster. Although in the fol- 
lowing we model the pulsar wind by a flow with the non- 
rclativistic speed of 10 4 kms~ , the crossing time for such 
a slow wind is only ~ 1 d. Therefore, the starting time 
docs not matter as long as it is at least several tens of 
days before periastron. 

In the tidal simulation, we emulated the Shakura- 
Sunyaev viscosity, adopting variable aspn over space and 
time and fixing /?sph to 0. In the wind simulation, how- 
ever, this method would fail, allowing the particle inter- 
penetration at strong shocks. Therefore, we adopted the 
standard values of the artificial viscosity parameters in the 
wind simulation, i.e., «sph = 1 and /?sph = 2. In the wind 
simulations, the energy equation is also different from the 
tidal simulation. In the latter, the gas was assumed to 
be isothermal, in order to emulate the temperature dis- 
tribution in the Be disk. However, the purpose of the 
wind simulation is to study the structure of the wind-wind 
and wind-disk collisions. Therefore, in the wind simula- 
tion, we take account of optically thin radiative cooling. 
Numerical implementation of radiative cooling is done by 
adopting Townsend (2009) 's Exact Integration Scheme 
for radiative cooling, with the cooling function generated 
with CLOUDY 90.01 for an optically thin plasma with 
solar abundances (Fcrland 1996). 

With a large binary separation, the winds are likely to 
collide after the Be wind reaches a terminal speed. Thus, 
for simplicity, we assume that the winds coast without any 
net external forces, assuming in effect that gravitational 
forces are either negligible (i.e. for the pulsar wind) or 
are canceled by radiative driving terms (i.e. for the Be 
wind). The relativistic pulsar wind is emulated by a non- 
rclativistic 10 4 kms _1 wind with the adjusted mass-loss 
rate so as to give the same momentum flux as a relativistic 
flow with the assumed energy, as in Romero et al. (2007). 

In the following, for simplicity, we assume that all 
the spin down energy .Epsr, = 8.2 x 10 35 erg s _1 goes 
to the kinetic energy of a spherically symmetric pulsar 
wind. We also assume the Be wind to be spherically 
symmetric. Given a recent high-resolution spectroscopic 
study of LS 2883 suggesting a significantly earlier spec- 
tral type than the conventionally used spectral type B2V 
(Ncgueruela et al. 2011), we take two different values of 
mass-loss rates, 10 _9 M Q yr _1 for the conventional spec- 
tral type B2V and lO _8 M0yr _1 for an earlier spectral 
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Fig. 3. Snapshots over a month around periastron passage for po = 10 _11 gcm~ 3 , where po is the base density of the Be disk. The 
upper panels show the column density along the z-axis (binary orbital axis) in cgs units, while the lower panels the column density 
along the y-axis (minor axis of the binary orbit). Annotated in each panel are the phase and time measured from periastron passage 
and the number of SPH particles. 



type. The mass injection rate to the Be disk is fixed to 
3.5 x lO^Mgyr" 1 . 

3.2. Interaction between the pulsar wind and the Be disk 
and wind 

Once winds are turned on, the faster pulsar wind soon 
fills up the whole simulation volume of r < a. Then, 
the slower Be wind pushes the pulsar wind back to the 
radii where the ram pressures of both winds are balanced. 
Without the Be disk, the shape of the interaction surface 
would be determined by the wind-wind collision. A key 
parameter would, then, be the ratio of wind momentum 
fluxes, given by 

-Epsr 

V= Tr v ' () 

where V wln d and M W ind are the velocity and mass loss rate 
of the Be wind, respectively, and £psr is the power of the 
pulsar wind. Taking £psr = 8.2 x 10 35 erg s _1 for the 
spherically symmetric pulsar wind and V w i n d = 10 3 kms - 
as the velocity of the Be wind, we have rj ~ 4.3 for the 
"weak" wind case (M win d = 10 _9 M @ yr _:l ) and ?; ~ 0.43 
for the "strong" wind case (M W ind = 10~ 8 M©yr _1 ). Note 
that the pulsar wind dominates the Be wind in the former 
case, while the Be wind is about twice stronger than the 
pulsar wind in the latter case. Thus, having an accurate 
spectral type of the Be star is important to construct any 
satisfactory model for this VHE gamma-ray binary. 

In a simple 2D model that ignores orbital motion, the 
ram pressure balance then implies that, for a binary sep- 
aration D, the interface should be located at a distance 

, _ D ( Q.60D (strong wind) , > 

1 + 07 ~ \ °- 33D ( wcak willd ) 
from the Be star (Stevens et al. 1992; Canto et al. 1996). 



In adiabatic shocks, the interaction surface approaches a 
cone with the half opening angle, 

.-»(t--„r~{i°r g"? m 

(Gayley 2009), measured from the Be star. Thus, the 
shock is wrapped around the Be star in the weaker wind 
case, while it is around the pulsar in the stronger wind 
case. 

In PSR B1259-63/LS 2883, the presence of the Be disk 
adds extra complications in the interaction with the pulsar 
wind. Figure 4 shows how the interaction occurs around 
periastron in the weak (upper figure) and strong (lower fig- 
ure) wind simulations. In each figure, the upper and lower 
panels respectively present the column density along the 
z-axis (binary orbital axis) along the y-axis (minor axis 
or the binary orbit). The change in the Be disk structure 
from that in the tidal simulation (Fig. 3) is drastic. The 
pulsar wind now strips off an outer part of the Be disk on 
the side facing the pulsar, truncating it at a radius where 
the gas pressure of the disk is roughly comparable to the 
ram pressure of the pulsar wind. Since the pulsar wind 
has a velocity component tangential to the line connect- 
ing the Be star and the pulsar, the material in the disk 
outside this truncation radius is pushed sideways, forming 
a warped filament in the prograde direction. This feature 
is more remarkable in the weak wind simulation than in 
the strong wind simulation, because, unlike in the latter 
simulation, the Be wind in the former simulation is too 
weak to shield the disk from the pulsar wind. 

The interaction between the pulsar wind and the Be- 
star circumstcllar environment is more clearly seen in 
Fig. 5, in particular in the upper panels, where the vol- 
ume density in the binary orbital plane (upper panels) 
and along the disk mid-plane (lower panels) is shown at 
the same phases as in Figs. 3 and 4. Comparing Fig. 5(a) 
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(a) M wind = lO- 9 M yr- 1 
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(b) M wind = lO- 8 M yr- 1 
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Fig. 4. Snapshots over a month around pcriastron passage from simulations in which both the pulsar wind and the Be wind are taken 
into account: (a) A'/ W j n< j = 1O -9 M0 yr — 1 and (b) M w i n d = 10 — 8 Af©yr . The mass-loss rate via the Be disk is 3.5 X 10 Af©yr , 
where M w i n( j is the mass-loss rate through the Be wind. In each panel, Ni, N2, and JV3 annotated at the lower left corner are the 
numbers of particles in the Be wind, the pulsar wind, and the Be disk, respectively. 



for the weak wind case with Fig. 5(b) for the strong wind 
case, we note that the effect of the pulsar wind on the Be 
disk depends on the relative strength of the Be wind. In 
the case of weak Be wind, the Be disk is not only truncated 
but also strongly deformed by the pulsar's ram pressure 
around pcriastron. In contrast, little deformation is seen 
in the case of strong Be wind, where the Be wind shields 
the Be disk from the pulsar wind. 

In both simulations, the truncation of the Be disk is 
so efficient that the disk has a sharp density drop at the 
outer radius. At phases of closest encounter of the pulsar 
with the disk outer radius, the density at the outer radius 
is enhanced by the shock by a factor of several, as seen in 
the lower panels of Fig. 5. This sharp truncation of the 



Be disk is likely a characteristic to be seen only in systems 
consisting of a Be star and an object with a strong wind. 
Although Be disks in binaries are, in general, tidally trun- 
cated, the density decrease beyond the disk outer radius 
is much more gradual in the tidal truncation. The density 
drop near the central star is due to ablation by the Be 
wind. 

Despite these features related to the Be disk, the loca- 
tion and global shape of the interaction surface between 
the pulsar wind and the circumstcllar environment of the 
Be star, i.e., the Be disk and wind, is consistent with that 
expected in the above analysis using a simple 2D model 
without the Be disk. We can see that in Fig. 5, where the 
density on the binary orbital plane is shown at the same 
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phases as in Figs. 3 and 4, the apex location and coni- 
cal shape of the interaction surface in the simulation, in 
general, agree with those given by equations (2) and (3). 

4. Summary and Discussion 

In this paper, we have numerically investigated the hy- 
drodynamic interaction of the circumstellar environment 
of the Be star with the pulsar in the TeV binary PSR 
B1259-63/LS 2883. We have broken the computation 



into two separate, but linked parts. The first part con- 
siders only the gravitational interaction with the pulsar, 
ignoring the effects of both the pulsar and Be star winds. 
One purpose of this part of simulation is to see how the 
Be disk evolves by the effect of viscosity and modulates by 
the tidal interaction with the pulsar. The other purpose 
is to set up the initial configuration for the second part 
of simulation, which takes account of both pulsar and Be 
winds. 

From the first part, we confirmed the viscous decretion 
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disk scenario that the material ejected from the stellar 
equatorial surface drifts outward by the effect of viscosity, 
and forms a Keplerian disk. Note that the radial velocity 
caused by viscous diffusion is very subsonic, so it takes 
more than a decade for this very wide system to have a 
fully developed Be disk. Without the pulsar wind, the Be 
disk is not truncated in such a highly eccentric system, 
unlike tidally truncated Be disks in binaries with lower 
eccentricities. Then, the pulsar passes through the disk 
twice an orbit, causing the tidal stream between the disk 
and the pulsar. 

The second part focuses on the interaction between the 
pulsar wind and the circumstcllar material of the Be star 
around periastron. It shows that once winds from the 
pulsar and the Be star are switched on, their effect soon 
becomes apparent in the structure of the Be disk. The 
pulsar truncates the Be disk on the side facing the pulsar, 
sweeping up an outer part to a dense filament. The result- 
ing structure of the disk is strongly asymmetric and phase 
dependent. The size of the truncated disk is so small that 
the pulsar never passes through the disk. Despite these 
features, the location and global shape of the interaction 
surface between the pulsar wind and the Be wind is well 
described by a simple 2D model. 

We performed numerical simulations of Newtonian hy- 
drodynamics without magnetic fields, although both of 
relativity and magnetic fields can play important roles on 
the dynamics of the system. A highly relativistic pulsar 
wind was emulated by a Newtonian wind by equating its 
momentum flux (thrust). Some studies of special relativis- 
tic hydrodynamic simulations show that highly supersonic 
flows display extended cocoons with high pressure (Marti 
et al. 1997), which correspond to the shocked pulsar 
winds in this study. The cocoons (or pulsar winds) have 
smoother local structure in relativistic simulations than in 
non-relativistic simulations, because the inertia of a rela- 
tivistic flow is smaller than that of a non-relativistic flow 
with the same thrust (Rosen et al. 1999). The global 
structure of the shocked pulsar wind, however, will not 
considerably be changed, since we have tuned the wind 
momentum flux in the simulation identical to that of rel- 
ativistic wind of PSR B1259-63. The overall structure of 
the shock is basically determined by the balance of mo- 
mentum fluxes of the flows. It is unlikely as well that 
the magnetic fields strongly affects the overall interaction 
feature, given that the kinetic energy is, in general, much 
greater than the magnetic field energy in pulsar winds 
near termination shocks (e.g. Kennel & Coroniti 1984). 
Introduction of magnetic fields into the simulation, how- 
ever, may be needed to study the local structure in detail 
(McKinney 2006; Kommisarov et al. 2009; Nagataki 
2009; Nagataki 2010). Therefore, it will be important to 
extend our code to a relativistic regime (e.g. Monaghan & 
Price 2001; Ryu et al. 2006; Rosswog 2010) with mag- 
netic fields (e.g. Price & Monaghan 2005) in the future 
so that we can self-consistcntly study these effects. 

It is frequently discussed that the particle acceleration 
happens at the shock in the pulsar wind side and the 
accelerated, non-thermal electron-positron pairs emit ra- 



dio to X-ray photons through synchrotron radiation while 
VHE (GeV-TeV) gamma-rays are produced by inverse- 
Compton scatterings of soft photons from the Be star 
(Tavani & Arons 1997; Sicrpowska-Bartosik & Bednarek 
2008; Takata & Taam 2009; Naito et al. 2010). It is 
also claimed that proton acceleration at the shock in the 
Be star side and VHE gamma-rays from pion decays that 
are produced through proton-proton interactions may be 
relevant to explain VHE gamma-rays (Kawachi et al. 
2004; Chernyakova et al. 2006). From the point of view 
of the time-correlation among radio, X-rays, and VHE 
gamma-rays (Chernyakova et al. 2006), hadronic scenario 
may be favored. This is because the anti-correlation of 
fluxes between synchrotron emission and inverse-Compton 
scatterings is expected for the leptonic scenario, which 
seems to contradict the observations (Tavani & Arons 
1997). To explain the observed correlation between X-ray 
and VHE emissions, therefore, Takata & Taam (2009) 
invoked the scenario that the observed emission in differ- 
ent orbital phases emanate from different magnetic field 
lines, giving rise to the pulsar wind parameters chang- 
ing with the orbital phase or that the power law index 
of the accelerated particles varies with the orbital phase. 
However, since these discussions depend on one-zone mod- 
els adopted, the anti-correlation may not appear even for 
the leptonic scenario when we investigate a more realistic 
situation using the results from our simulations. For ex- 
ample, the previous leptonic models have assumed that 
the pulsar is confined by the shock (Tavani & Arons 
1997; Takata & Taam 2009), whereas the result from 
the present simulations suggests that the shock geometry 
is sensitive to the mass loss rate from the Be star and to 
the orbital phase, as Figure 5 shows. It is thus essential to 
calculate the high energy emission processes with a more 
realistic situation. We are planning to investigate it as 
our next step. 

While a pulsar wind + stellar wind interaction model 
has also been proposed for other high-mass binaries with 
GeV-TeV gamma-ray emission (e.g., LS I +61 303 and 
LS 5039), PSR B1259-63 is so far the only such system 
with direct detection of pulsed radio emission. Johnston et 
al. (1996) showed that this radio emission disappears for a 
period of about 5 weeks centered on the periastron epoch 
of 1994 January 9. They attributed this as "most likely 
due to a combination of free-free absorption and severe 
pulse scattering in the Be-star disc" . But note that GHz 
radio emission is expected also to undergo strong free-free 
absorption by the much lower-density Be-star wind, since 
for typical wind parameters, the associated free-free opti- 
cal depth can be of several thousand at orbital distances 
around an AU [see, e.g., equation (4) of Torres (2010)]. As 
such, the detection of radio pulses in PSR B1259— 63 may 
actually be due to an observer perspective that is nearly 
aligned to the major axis of the binary orbit, such that 
during most of the orbit around apastron, the observed 
radio emission propagates mostly through the low-density 
shock cone surrounding the relativistic pair wind. In this 
scenario, the disappearance of radio emission near perias- 
tron could actually be attributed to the rapid rotation of 
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the shock cone away from the viewer's line of sight, caus- 
ing strong free-free attenuation of the radio emission by 
the Be-star wind and its associated shock cone. A similar 
viewer perspective has been used to explain the several- 
week-long attenuation of X-ray emission in the colliding 
wind binary 77 Carinae (Okazaki et al. 2008). Our future 
work will further explore this shock-cone scenario as an 
alternative to the standard Be-disk attenuation model for 
the periastron disappearance of pulsar radio emission. 

Observations of TeV binaries in high energy bands pro- 
vide important information on the type of interaction and 
the emission mechanisms. For TeV binaries with Be stars, 
however, the high energy bands are not the sole window to 
probe these systems. As mentioned in section 2, the op- 
tical emission lines and infrared excess arise from the Be 
disk, while the Be wind emits UV radiation. Observations 
in these bands are thus suitable to study the effects of the 
pulsar wind on the circumstellar environment of the Be 
star. Particularly, spectroscopic and polarimctric observa- 
tions in the optical and infrared can provide valuable infor- 
mation on the structure and dynamics of the disturbed Be 
disk, which gives clues for the interaction between the pul- 
sar wind and the Be disk. UV observations, on the other 
hand, are adequate for studying the strongly asymmetric 
Be wind around periastron (see Fig. 5 for the asymmet- 
ric structure of the Be wind). For better understanding 
of TeV binaries with Be stars, including PSR B1259— 63, 
multi-band observations are highly desirable. 
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